/*******************************************************************************
 * Copyright 2009 DCSpectrometer - http://code.google.com/p/dcspectrometer 
 *  
 * Licensed under the Apache License, Version 2.0 (the "License"); 
 * you may not use this file except in compliance with the License. 
 * You may obtain a copy of the License at 
 *  
 *     http://www.apache.org/licenses/LICENSE-2.0 
 *     
 * Unless required by applicable law or agreed to in writing, software 
 * distributed under the License is distributed on an "AS IS" BASIS, 
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 
 * See the License for the specific language governing permissions and 
 * limitations under the License. 
 *******************************************************************************/
package com.dcspectrometer.Analysis;

/**
 * This class calculates Cumulative Distribution Function for Normally Distributed Variables
 * Uses a 5-point estimation formula.
 */
public final class Norm5PointCDF {
	private static final double coefficient = (1/(Math.sqrt(2*Math.PI)));
	/**
	 * Calculates probability that normal random variable is less than x
	 * @param x - input value
	 * @return - Probability that normal random variable is less than x
	 */	
	public static double getCDF(double x) {
		final double a1 = 0.3193;
		final double a2 = -0.35656;
		final double a3 = 1.78147;
		final double a4 = -1.82125;
		final double a5 = 1.33027;
		final double k1 = 1/(1+0.2316419*x);
		final double k2 = k1*k1;
		final double k3 = k2*k1;
		final double k4 = k3*k1;
		final double k5 = k4*k1;
		return 1 - (coefficient*Math.exp(x*x/(-2)) * (a1*k1 + a2*k2 + a3*k3 + a4*k4 + a5*k5));
	}
	
	/**
	 * Computes a derivative for a normal random variable CDF for a given x 
	 * 
	 * @param x - input value
	 * @return - derivative at value
	 */
	public static double getCDFPrime(double x) {
		final double a1 = 0.3193;
		final double a2 = -0.35656;
		final double a3 = 1.78147;
		final double a4 = -1.82125;
		final double a5 = 1.33027;
		final double y = (1+0.2316419*x);
		final double y2 = y*y;
		final double k1 = 1/(1+0.2316419*x);		
		final double k1prime = -1 * 0.2316419/(y2);
		final double k2 = k1*k1;
		final double k2prime = 2 * k1prime/(y);
		final double k3 = k2*k1;
		final double k3prime = 1.5 * k2prime/(y);
		final double k4 = k3*k1;
		final double k4prime = 2 * k2prime/(y2);
		final double k5 = k4*k1;
		final double k5prime = 1.25 * k4prime/(y);
		final double v = (a1*k1 + a2*k2 + a3*k3 + a4*k4 + a5*k5);
		final double vprime = (a1*k1prime + a2*k2prime + a3*k3prime + a4*k4prime + a5*k5prime);
		final double u = coefficient*Math.exp(x*x/(-2));
		final double uprime = -coefficient*x*Math.exp(x*x/(-2));
		
		return - (v * uprime + u * vprime);
	}
}